# convert felm result into a data frame with rhs | coef | se | xmin | xmax | xmid | ci.l | ci.h,
# replacing Inf and -Inf with in xmin and xmax on the edge bins with the graphically spaced bound
## felm.est: result from felm() call (lm will probably work too) with RHS factor variable like cuttmax or bins
## plotvar: string of variable to plot: e.g. 'cuttemp'
## breaks: bin size
## omit: omitted bin

bin.plot <- function(felm.est, plotvar, breaks, omit, panel="panel", ylabel="y-label", xlabel="x=label",
                     y.axis.title = element_text(size=30, angle=90, margin = margin(t = 0, r = 10, b = 0, l = 0), face="bold"),
                     ylimit = c(min(dt$ci.l) - (max(dt$ci.h) - min(dt$ci.l))/10,max(dt$ci.h)),
                     y.axis.text = element_text(size=25, angle=0, face="bold"),
                     y.axis.line = element_line(),
                     y.axis.ticks = element_line(),
                     errorline="#8da0cb", errorfill="transparent", linecolor="#fc8d62", pointfill="red4", pointcolor="#8da0cb") {

    dt <- as.data.frame(summary(felm.est)$coefficients[, 1:2]) # extract from felm summary (use df to keep coefficient names)
    dt <- dt[grepl(plotvar, rownames(dt)), ] # extract variable name (e.g., tmax.cut, ppt.cut, ...), limit to matches
    dt$rhs <- str_split_fixed(rownames(dt), "   ", n=2)[,1]
    dt <- as.data.table(dt)
    dt[, rhs := gsub(paste0(plotvar), "", rhs)]
    dt[, rhs := gsub("\\(|]", "", rhs)]
    dt[, rhs := gsub("neg", "-", rhs)]
    dt[, rhs := gsub("pos", "", rhs)]
    dt[grepl(",", rhs), rhs := tstrsplit(rhs, ',')[2]] # if a factor variable, grab second column (upper limit)
    dt[, rhs := as.numeric(rhs)]
    dt <- dt[!is.na(rhs)]
    dt[, xmin := rhs - breaks]
    dt[, xmax := rhs]
    dt[, rhs := NULL]
    setnames(dt, c("coef", "se", "xmin", "xmax"))
    dt[, xmin.lead := data.table::shift(xmin, type="lead")]
    dt[, xmax.lag := data.table::shift(xmax, type="lag")]
    dt[xmin==-Inf, xmin := xmin.lead - breaks]
    dt[xmax==-Inf, xmax := xmin.lead]
    dt[xmin==Inf, xmin := xmax.lag]
    dt[xmax==Inf, xmax := xmin + breaks]
    dt[, c('xmin.lead','xmax.lag') := NULL]
    # identify omitted category, add 0.
    dt <- rbind(dt, data.table(coef=0, se=0, xmin=omit[1], xmax=omit[2]))
    dt[, ':='(xmid=(xmin+xmax)/2, ci.l=coef - se*1.96, ci.h=coef + se*1.96)] # add some stuff for plotting
    dt[, range := paste0(xmax)] # create range variable for x labels
    dt[xmax == min(xmax), range := paste0("<=",min(xmax))] # set bottom of range
    dt[xmin == max(xmin), range := paste0(">",min(xmin))] # set top of range

    plot <- ggplot() +
        geom_hline(aes(yintercept=0), linetype=2, color='gray50') +
        geom_ribbon(data=dt, aes(x=xmid, ymin=ci.l, ymax=ci.h), fill=errorfill, alpha=0.1) +
        geom_line(data=dt, aes(x=xmid, y=ci.l), linetype="dotted", color=errorline, size=1.5, alpha=1) +
        geom_line(data=dt, aes(x=xmid, y=ci.h), linetype="dotted", color=errorline, size=1.5, alpha=1) +
        geom_line(data=dt, aes(x=xmid, y=coef), color=linecolor, size=1.75, alpha=1) +
        geom_point(data=dt, aes(x=xmid, y=coef), size=3, fill=pointfill, color='black', pch=21, alpha=1) +
        annotate("text", label=panel, x = min(dt$xmid) + (max(dt$xmid) - min(dt$xmid))/10, y = min(dt$ci.l) - (max(dt$ci.h) - min(dt$ci.l))/10, size=9) + # plot panel annotation
        ylab(ylabel) +
        xlab(xlabel) +
        coord_cartesian(xlim=c(min(dt$xmid), max(dt$xmid)), ylim=ylimit, expand=TRUE) +
        scale_x_continuous(breaks = dt$xmid, labels = dt$range) +
        theme(plot.title = element_text(face="bold", size=40),
              axis.title.x = element_text(size=30, angle=0, face="bold"),
              axis.title.y = y.axis.title,
              axis.text.y = y.axis.text,
              axis.text.x = element_text(size=25, face="bold"),
              panel.grid.major = element_blank(),
              panel.grid.minor = element_blank(),
              panel.background = element_blank(),
              axis.line.x = element_line(),
              axis.line.y = y.axis.line,
              axis.ticks.y = y.axis.ticks,
              legend.title=element_blank()
        )
    return(plot)
}
